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We describe a time-resolved monitoring technique for heterogeneous media. Our approach is based on the 
spatial variations of the cross-coherence of coda waveforms acquired at fixed positions but at different dates. 
To locate and characterize a weak change that occurred between successive acquisitions, we use a maximum 
likelihood approach combined with a diffusive propagation model. We illustrate this technique, called LOCAD- 
IFF, with numerical simulations. In several illustrative examples, we show that the change can be located with 
a precision of a few wavelengths and its effective scattering cross-section can be retrieved. The precision of the 
method depending on the number of source receiver pairs, time window in the coda, and errors in the propaga- 
tion model is investigated. Limits of applications of the technique to real-world experiments are discussed. 



I. INTRODUCTION 

Elastic and acoustic waves constitute one of the primary 
tools to detect and locate temporal changes in natural or man- 
made structures. If the waves do not interact with any other 
obstacle than the target, conventional imaging techniques are 
usually based on geometrical considerations. A controlled 
pulse emitted into the medium is scattered by the target and 
the echos are recorded with a receiver. These techniques can 
be improved using several sources and detectors, and extended 
to locating several targets at the same time. As long as the typ- 
ical propagation time in the medium is much smaller than the 
scattering mean free time, i.e. the average time between two 
scattering events, we are in the single scattering regime. In 
this case, the resolution for detecting and locating a change is 
limited by the Fresnel zone yXL, with L the typical propa- 
gation distance in the medium. Applications in every day life 
abound: they cover high-stake fields like ultrasonic medical 
imaging, non-destructive testing, seismic exploration, radar 
aircraft location or sonar. 

This simple picture does not apply in heterogeneous me- 
dia such as polycrystals, concrete, or volcanoes. Imaging 
these materials in a non-destructive way is an important is- 
sue for miscellaneous applications like monitoring, ageing or 
damage assessment. In heterogeneous media, ray theory is 
not relevant because the scattering mean free time is much 
smaller than the typical record duration. A pulse emitted into 
the medium experiences numerous scattering events and the 
output signal recorded at large distance from the source dis- 
plays complex details that depend on the interactions between 
the wave and each of the scatterers. Beyond a distance called 
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transport mean free path £*, the memory of the initial direc- 
tion of propagation is lost. In this regime, the average energy 
distribution in the medium evolves as a diffusion process and 
it is relevant to describe wave propagation using probabilities. 

The problem of locating an isolated change in a multiple 
scattering sample has received some attention in the past, par- 
ticularly in optics. The space and time correlations of inten- 
sity in a speckle pattern probed by one or more receivers al- 
low one to observe the diffusion of scatterers iQly]. On one 
hand, diffusive wave spectroscopy yj] and its variants have be- 
come standard tools for investigating collective changes in the 
medium. On the other hand previous authors |4] have shown 
that a local change (the perturbation) within a collection of 
scatterers (the background) essentially acts as a dipole source 
of intensity. Intensity variations enable the detection and lo- 
cation of a crack from observations in transmission |5l|g]> or 
more generally to locate an object with known characteris- 
tics J3>[Ml- The weak sensitivity of the method has been il- 
lustrated by numerical studies [6]. Indeed, a large amount of 
ensemble or frequency averaging (typically 100 realizations) 
is required to distinguish the intensity fluctuation caused by 
the defect from the background speckle pattern. From a theo- 
retical point of view, the weak sensitivity can be traced back 
to the cancellation of diagrams that dominate the waveform 
decorrelation, a cancellation which is imposed by the optical 
theorem. This renders techniques based on intensity varia- 
tions almost inapplicable to solid media. These points will be 
further illustrated below. 

In acoustics, one can commonly record a large number of 
signals with perfect temporal and spatial resolution, which 
is advantageous compared to optics. A pulse emitted into a 
medium gives rise to long time records with a pronounced 
coda, a term which refers to the arrivals following the ballistic 
pulse. Several techniques use the coda to retrieve information 



on the evolution of the medium. In seismology, the monitor- 
ing of temporal changes in the crust was initiated in the mid- 
80's, using repeating small earthquakes on faults fl. Later on 
the method was applied to volcanoes and revealed temporal 
changes of velocity prior to eruptions [10]. The method was 
transposed to the laboratory and popularized under the terms 
diffuse acoustic wave spectroscopy (DAWS) lllil . or coda- 
wave interferometry (CWI) 11121 1 1 311 . In these approaches, 
changes of waveforms in the coda are interpreted in terms of 
travel time variations, a technique that is very sensitive 111411 
to detect weak changes, but gives little information concern- 
ing the location of the change. To first order, global velocity 
changes in the medium result in a stretching of the waveforms 
Ill0l[l4l - [l6ll but the interpretation of a local change in terms of 
travel time fluctuation remains problematic. Recently DAWS 
has been used in damage monitoring [ 17, 18] but a large range 
of other applications are possible 11911 . For a broad review of 
applications of CWI in geophysics, we refer to [20]. Also 
based on the concept of correlation, techniques have been de- 
veloped to recover the Green's function in an open medium 
based on the cross correlation of noise signals J2ll425ll . These 
noise-based Green's functions can in turn be used in a pas- 
sive image interferometry technique with applications in vol- 
canology and fault monitoring Ilia, uol un\ . Recently, Aubry 
and Derode [28] proposed an alternative technique based on 
the singular value decomposition of the propagator, but this 
technique is limited to a sufficiently strong extra scatterer and 
is not sensitive to weak perturbations. 

In this article, we report on a different approach to locate 
a small isolated change. Our LOCADIFF technique uses the 
correlations between time windows in the late coda for sev- 
eral pairs of sources and receivers. A numerical model of the 
medium is then used to compute the most likely position of the 
weak change, in terms of probability. We start our descrip- 
tion of the work by observing the correlation loss of signals 
induced by the weak change in a finite difference numerical 
simulation (Section HJi. Using the theory of multiple scatter- 
ing [8], we derive an expression of the decorrelation induced 
by a weak change in Section [III] We then present the inver- 
sion technique, based on the maximum likelihood principle in 
Section HV Al We discuss the accuracy of the technique and 
possible improvements in Section lVTl 



A. Numerical simulations of wave propagation 

As a first investigation, we perform 2D numerical exper- 
iments of acoustic wave propagation in heterogeneous open 
media [i29il . Using a finite difference centered scheme, we 
solve the wave equation with absorbing boundary conditions; 
the dimension of the simulation grid is 50Ao X 50Ao, with a 
spatial discretization step Ao/30, where Ao is the central wave- 
length. Synthetic data are computed on a linear array of 9 re- 
ceivers located at the center of the medium and 10 sources are 
randomly distributed over the grid. Sources and receivers are 
kept fixed throughout the experiments (see Fig.Q]). To mimic 
a multiple scattering medium, 800 empty cavities of diameter 
Ao/3 are randomly distributed over the grid. In the frequency 
band of interest the average scattering cross-section was nu- 
merically estimated as £ = 1.6Ao, along with the transport 
cross section E* = l.lArj. Table U summarizes the physi- 
cal properties of the simulated medium, including the number 
of scatterers (with density n), the transport mean free path 
£* = — iy, the diffusion constant D = £=- and the Thou- 

less time td — J7J, where R 2 is the mean squared distance 
between sources and receivers. Note that these quantities are 
evaluated under the "independent scattering approximation", 
which assumes that the waves never visit the same scatterer 
twice. 
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FIG. 1. Distribution of sources, receivers, scatterers in the numerical 
simulation. The scatterer to remove is surrounded in gray. 



II. OBSERVATIONS OF CORRELATION LOSS AFTER A 
WEAK CHANGE 



It is already known that a weak change can be detected in 
a scattering medium because it slightly modifies the coda of 
the Green's functions. The amount of modification is usually 
quantified by measuring the cross-correlation between wave- 
forms recorded at different times [9]. We illustrate the signal 
processing with the aid of finite-difference simulations of the 
wave equation in a medium containing a large number of iden- 
tical scatterers. 



Parameter 


notation 


value 


Number of scatterers 




800 


Transport mean free path 


t 


2.8 Ao 




kt 


18 


Diffusion constant 


D 


lAXl/T 


Thouless time 


TD 


68 To 


Coda decay time (leakage) 


T a 


240 T 



TABLE I. Physical parameters of the simulations in normalized units. 



The signal e(t) emitted by each source is a pulse with cen- 
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of duration 2T centered on t using the formula: 
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FIG. 2. Normalized waveforms h(t) obtained at receiver 1, for 
source number 1 (top) and 3 (bottom). Inset: zoom into the coda. 
The black solid (resp gray broken) line corresponds to the record ac- 
quired before (resp. after) the change. For illustration purpose, the 
value of To was set to 1 /^s. 



tral frequency /o and a Gaussian envelope (100% bandwitdh 
at -6dB). Using source i we record with receiver j the sig- 
nal hij (t) during 300 oscillations of period To. Typical wave- 
forms hij (t) are plotted in Fig. [2] The long tail of the record 
in Fig.|2]corresponds to arrival of partial waves that have been 
scattered several times. Notice that the ballistic arrival is not 
distinguishable in the waveforms of Figure [2] Long coda and 
lack of ballistic arrival constitute evidences that we are in a 
strongly scattering regime, in agreement with our estimates of 
the transport mean free path. During the first run of the simu- 
lation, 10 x 9 impulse responses hij are recorded and stored. 
On a second run, one scatterer is removed and another set of 
impulse responses h'^ is evaluated. Both hij(t) and h^At) 
display long codas lasting a large number of ballistic times. 
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The typical width of the time window T is of the order of 
5Tq. Experimentally, enlarging T partly eliminates the effect 
of noise and reduces the fluctuations of the correlation coef- 
ficient. Nonetheless, using a large value for T results in con- 
sidering simultaneously paths with very different lengths. We 
address this important point in Section lTllBI 



C. Spatial dependence of the decorrelation 

In Figure |2j it is noticeable that the differences between 
hi,i(t) and h[ x (t) (top), are much smaller than the dif- 
ferences between /i3,i(t) and h' 3 1 (t) (bottom), even in the 
late coda. The decorrelations computed over the interval 
[210 T , 220 T ] are A'i,! (215 T ) = 5% and #3,1 (215 T ) = 
27%, respectively. Consequently, the amount of decorrela- 
tion depends on the positions of the source and receiver with 
respect to the local change, a property which holds even in 
very late time windows in the coda. For a given configuration 
of source -receiver pairs, we obtain a set of observed decor- 
relations, which are characteristic of the relative locations 
of the source, receiver and change in the multiple scattering 
medium. We will now demonstrate the possibility to locate the 
change and estimate its cross-section from the knowledge of 
the source and receiver positions and the corresponding decor- 
relation coefficients. To do so, we develop a theoretical model 
to predict the decorrelation coefficient of waves induced by 
the addition of a change in a heterogeneous medium, in the 
diffusive regime. We recall in the next section the necessary 
elements from multiple scattering theory. 



B. Detection of an isolated change 



III. WAVE SCATTERING THEORY 



We assume that the medium can be represented as a ma- 
trix with embedded inclusions. Only the scalar case is con- 
sidered here. The scattering properties of an inclusion will 
be described by its T matrix, defined in operator notation as 
MM: 



G\ — Go + GqTGq 



The details of the complex waveforms shown in Figure [2] 
are highly sensitive to the positions of the the scatterers. Each 
waveform can be understood as a fingerprint of the medium. 
As our goal is to detect a single scatterer's change, we need 
to exploit the information contained in both the amplitude and 
phase of the signals. A comparison of the records hij (t) and 
h'^(t) reveals that for short coda times up to 100 T « 30r*, 
no difference is visible in the signals. We observe small dif- 
ferences between the waveforms at later times that are solely 
due to the change in the medium. Figure [2] (bottom) shows 
an example of such differences for hz,iit) and /ig x (t). The 
observed decorrelation is too large to be attributed to numer- 
ical noise, there is thus evidence that the coda waveforms are 
sensitive to the removal of only one scatterer. 

The differences between the waveforms hij and h' t , are 
quantified by the decorrelation, or correlation loss, between 
hij and h[ _■• The decorrelation is computed in a time window where a is the scatterer cross-section. 



(2) 



where Go is the retarded free space Green function and G\ is 
the Green function in the presence of the scatterer. For a non 
absorbing scatterer, energy conservation implies the following 
optical theorem: 
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FIG. 3. The Bethe-Salpeter equation defining the so-called ladder 
operator L. The solid line and dashed line represent the retarded 
ensemble averaged Green's function and its complex conjugate, re- 
spectively. The dotted line connecting the two vertices indicates that 
they represent the same scatterer. 



A. Correlations between two slightly different media 

We want to predict the decorrelation of waveforms in a 
medium where a small change occurs. Although we will em- 
ploy a statistical approach based on ensemble averages, in 
general we have access to only one realization of the random 
process. Therefore we introduce the following estimator of 
the cross-correlation function based on the observation of a 
single coda: 



1 



t+T 



r(«, r) = — J i?(t' + r/2)^(t' - r/2)dt', (4) 

where ip is the scalar field. The superscript 2 refers to the 
medium in presence of an extra defect while the superscript 1 
refers to the medium without it. We have introduced an analog 
of the Wigner function which is most convenient to analyze 
non-stationary signals. The empirical cross-correlation can be 
decomposed over internal and external frequencies u and ft, 
respectively: 



r(t, r) = 



1 



dfi 



dwr(fi, u>) exp[-i(fM+cj-r)] 
(5) 



where the frequency-domain cross-correlation reads: 

f(0, u) = sinc(OT) ip 2 (uj + Q/2)^{u - fi/2)*. (6) 



We see that the width of the time window, IT, has a minor 
effect only. Equation (|6]l shows that we have to compute the 
quantity 



(G 2 (lu + Q/2)G 1 (uj - Q/2)*) 



(7) 



In equation Q, G is the retarded Green's function. We will 
denote by 7o the T-matrix of the additional defect which is 
assumed to appear at the position Xrj. In diagrammatic nota- 
tions, such as the one employed in Figure [3] , the T matrices 
are represented by crosses. The transport of energy in the scat- 
tering medium is described by the ladder operator L, which is 
defined by the diagrammatic self-consistent equation shown 

in FigureESlfiol 13211. 

We use the field-field correlation function in the coda 



r(w, Q, s, xrj, r) = / dri / dr 2 P (w, ft, s, n) 

x L e (ui, ft, ri, x , r 2 )P (w, fij r 2 , r). (8) 



Quantities labelled with " are implicitely evaluated at inner 
frequency u) and outer frequency Q. The ladder propaga- 
tor L e describes the transport of correlations in a sequence 
of scattering events in the medium with an extra scatterer. 
Po(s, ri) and -Po( r 2, r) describe the ballistic propagation 
from the source to the first scattering event, and from the last 
scattering event to the detector, respectively: 



Po(ri, r 2 ) 



-R/l 



(47Ti?) 5 



jnn/c 



(9) 



where R = |r 2 — i"i| and c = d^k^uj) is the group veloc- 
ity at the frequency uj. The ladder propagator with the extra- 
scatterer L e is related to the ladder propagator without the 
extra-scatterer L as follows [4]: 



x 



L, 



L 



FIG 4. The diagram of the ladder operator with an extra scat- 
terer. The extra scatterer is sandwiched between two ladder oper- 
ators. Note that we have neglected the possibility that ensemble av- 
eraged Green's functions connect the extra scatterer with the source 
and/or the receiver which assumes that it is located at least one mean 
free path away from all sources and receivers. 



L e (s, xrj, r) =L(s, r)+ 

/dri dr 2 L(s, ri)J(ri, x , r 2 )i(r 2 , r). (10) 

In Equation ([Tol l, also represented by the diagram depicted in 
Figure |4] the first term represents the scattering paths that do 
not see the defect, while the second term describes the paths 
that visit the defect once. As we are in a regime of weak inter- 
action between the field and the scatterer higher-order terms 
can be neglected. We define the operator J that connects the 
two ladders by 

J(r, x , r') = 

/dr x /dr 2 G(r, v 1 )T(r 1 , r 2 )G(r 2 , r')G(r, r')*, (11) 

where G denotes the ensemble averaged Green's function. In 
the mesoscopic regime, J is evaluated to lowest order of the 
small quantity l/(ko£) <€i 1 for a point scatterer: 



J(r, x , r') 



il 2 To 
87rfco 



^(r-xo)5 (3) (x -r'). (12) 



Inserting expression ( fT2l into equation ([8]l one obtains: 
f(s, x , r) = / dri / dr 2 P (s, ri)i(ri, r 2 )Po(""2, r) 



dri / dr 2 P (s. ri)£(ri, x 



il 2 % ~ r 



8irk 



L(xq, r 2 )P (r 2 , r), 



(13) 



where the first term is the diffuse intensity in the medium 
without extra scatterer and the second integral is an inter- 
ference term caused by the extra scatterer. In the slowly- 
varying envelope approximation, the integrals can be evalu- 
ated to give: 

£ 2 ~ £ 2 ~ i£ 2 % ~ 

T(s, x , r) = -^L(s, r) - ^L(s, Xo )— L(x , r). 

(14) 
In the diffusive regime, the propagator of the wave intensity 
in the medium filled with scatterers, Pd, is the solution of the 
following equation 



-m-DVDPJn, r 2 ) = *W(n 



r 2 



(15) 



where D is the diffusivity. The ladder L is related to Pd by 



Lir-i, r 2 ) 



~-Pd(ri, r 2 ). Using these notation the diffuse 



intensity for a unit point source satisfies: 

f(s, x , r) = —Pd(s, r) + -r-Pd(s, x ) 7 — ^F d (x , r). 
Air An 2ko 

(16) 
In order to obtain the correlation function in the time domain, 
we double invert the Fourier transform over the variables u 
and Q, We further assume that the signal has been filtered in a 
narrow frequency band Acj in which the scattering properties 
vary little. Upon integration over cj and application of the 
optical theorem ©, the correlation function for a unit point- 
source normalized by the bandwidth Aw reads: 

r(s, x , r, t) = P d (s, r, t) 

ca /"' 
-— / duP d {s, x , u)P d (x , r, t-u). (17) 
* Jo 

We have therefore obtained the theoretical decorrelation 
K(xo,t) = ^Q(s, x , r, t), where 



Q(s, x , r, t) 



J Q duP d {s, x , u)P d {xo, r,t-u) 
Pd(s, r, t) 



(18) 

The negative sign in ( fTTI i comes from the optical theorem (en- 
ergy conservation) and ensures that the cross-coherence is less 
than one. The derivation presented in this section does not de- 
pend on the form of Equation ( fT5l l. which means that solutions 
to a more accurate transport equation can be substituted to Pd- 
Note that for a resonant point scatterer, a can be substituted 
with X 2 /tt. 



B. Computation of the decorrelation formula 



If the medium is absorbing, the same issue arises. In media 
with a uniform absorption time kT 1 , the absorption affects the 
numerator of Q in dT8T > by a factor exp[— nu — n(t — u)] = 
cxp[— nt] and the denominator by a factor exp[— Kb], There- 
fore, uniform absorption effects cancel out in the normal- 
ized decorrelation function, which is a genuine advantage of 
the present technique. In the case where absorption is non- 
uniform, it will affect differently Pd(s, xo) and Prf(xo, r) and 
the observed decorrelation pattern may be partly ascribed to 
the spatial variations of absorption. Consider a medium with 
constant diffusivity D and absorption k. The solution of the 
diffusion equation (TT~5T > in an infinite d-dimensional medium 
is 



Pd(r ll r 2 , t) = 



(4irDt) d / [ 



■ exp 



— Kt — 



(r 2 - n) 2 



4Dt 



(19) 

In the case of a 3-D infinite medium, a usual Laplace trans- 
form calculation gives the exact result: 



Q(s, x, r, t) = 



1 



AttD 



1 1 



exp 



*R 2 



ADt 



(20) 

where we have introduce the notations s = ||s — xj|, r = 
|| r — x|| and R = ||s — r||. We observe that Q is a function 
with elliptic contour lines multiplied by simple poles located 
at s and r. Of course, if r = s, we recover the formula derived 
in Ref. [12] for an infinite medium. This formula is gener- 
ally not applicable under this form because the transducers are 
usually located at the surface of the system. However, if the 
boundary conditions are sufficiently simple, the formula (l20l > 
can be used as a building block to derive more complicated 
solutions, as shown in Section [Vl 

In formula ( ITTb we neglect two constraints. First, we as- 
sume that the change occurs at a minimum distance of the 
order of one mean free path from the source and the receiver. 
Second, we neglect the finite velocity of the wave, in other 
words, the contribution for times u,t — u < R/c in the inte- 
gral ( [T7l i should be removed. The contribution of short times 
u < R/c in ( ITTb is negligible as soon as ct ^> R > I*. The 
computation of the decorrelation coefficients Kij(t) must be 
done with T larger than a few oscillation periods of the wave. 
Using formula d20b , we can estimate the correction due to this 
averaging as a function of T/t. To do so, we compute the av- 
erage of (l2"0t on the interval [t — T, t + T] and divide by the 
value of Q at t. We obtain a curve of relative correction as a 
function of T/t which is independent of any other parameters 
and which is displayed on Figure [5] In most applications, the 
correction will be typically less than 10%. 



We observe that the decorrelation ( fT8l can be computed if 
the function Pd is known. In the general case where the diffu- 
sivity D depends on the position, the function Pd can only be 
numerically estimated, provided that the spatial dependence 
of D is known. In practice, the decorrelation coefficient can 
be reasonably rapidly computed if one assumes that the value 
of D is approximately uniform in the medium. We investigate 
the amount of variation for D in Section lTV El 



C. Intensity variations vs field correlations 

As recalled in the introduction, a number of investigations 
on the monitoring of complex media have focused on the de- 
tection of intensity variations induced by local changes of 
the scattering properties. We will show that in the diffusive 
regime, intensity variations are much less sensitive to local 
changes than field correlations. To do so, we calculate the 
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FIG. 5. Deviation of the average of Q on the time interval [t — T, t + 
T] with repect to Q(t). The correction grows more rapidly for large 
values of the argument in the exponential of l |20t , denoted by A in 
this Figure. 



perturbation of the ladder propagator induced by an extra- 
scatterer following the approach developed in reference [4]. 
In addition to the diagram depicted in Figure [4] two other di- 
agrams contribute to intensity variations: 1) a diagram with a 
single cross on the lower line and 2) a diagram with one cross 
on each line which are connected by dashed line. In the dif- 
fusive regime and for a non-absorbing defect, we obtain the 
intensity perturbation to lowest-order in q and q' in the form: 



SL[{s, x, r; t) 



£ 4 TT* 



+ oo 



an 



48-n- 2 J_ o0 2ir jjj R 3 xR 3 
L(q;fi)e iq -( r - x) (q- q')i(q'; 0)e iq '- (x - s ) e - im 



d 3 qd 3 q' 
(21) 



In the Fourier domain, the ladder propagator in the diffusion 
regime writes: 



L(q; n) 



Air 
(2ir) 3 P {q 2 £/3 - iltyc) ' 



(22) 



After integration over the wavenumbers q, q' and the fre- 
quency il, we obtain: 



5L' e (s, x, r; t) 



ac 



. e -I?/ADt x 



67r l/2 L) 3/2 < 3/2 

V s • (V r Q(s, x, r)) . (23) 



After calculation of the partial derivatives, we obtain the fol- 
lowing formula for the ladder perturbation induced by an extra 
scatterer: 



5L T e (s, x, r;i) = 



<rc 2 (r — x) • (s — x) 

487T 3 / 2 D 7 /2t5/2 r 2 s 2 

(r + s) 3 



r 3 + s 3 



rs 



2Dt 



-(r+s) 2 /4Dt 



(24) 



The intensity variation exhibits a characteristic pattern with 
positive and negative lobes, depending on the cosine of the 
angle between the source and receiver as seen from the ad- 
ditional scatterer. Even more important is the temporal de- 
pendence £~ 5 / 2 which is faster than the temporal decay of the 



ladder propagator between the source and receiver. As a con- 
sequence, the sensitivity to the local change decays like 1/t 
in the coda in sharp contrat to the field correlation which goes 
to a constant at large lapse time. This property justifies the 
popular use of field correlation functions to monitor temporal 
changes in evolving media. 



IV. THE INVERSION TECHNIQUE 

A. Maximum likelihood of the position 

In Section [HI] we have obtained an expression for the 
expected decorrelation as a function of the position of the 
change. The principle of the inversion technique is to compare 
a numerical model to experimental data. The change is found 
at the position where numerical and experimental decorrela- 
tion match best. The mismatch is measured by a standard 
least-squares cost function (x 2 )- The inversion technique con- 
sists in finding the position x and the cross-section a mini- 
mizing the function \ 2 - Such a technique is also often called 
a maximum likelihood method. Let us chose a set of sources 
Si (1 __ * __ n s) and a set of receivers Vj (1 < j < n r ), 
and call A^ the number of source-receiver pairs (in this case, 
N = n r n s ). There is no restriction on their positions, and in 
particular source and receiver can be located at the same posi- 
tion. We describe the technique at fixed time t in the coda. 




FIG. 6. The function Qij (x) in a infinite medium in dimension 3 
with constant diffusivity and absorption computed with formula d20t . 
The function is plotted in the plan containing the source, the receiver 
and the change. The value along the z axis (logarithmic) is the sensi- 
tivity to a change at the position in (a;, y). The two peaks correspond 
to the positions of the source and the receiver. The z-scale is loga- 
rithmic and arbitrary. 

The most restrictive assumption of our approach is that a 
single defect affects the experimental values of the decorrela- 
tion. The LOCADIFF inversion technique consists in retriev- 
ing the most likely position of this defect by introducing the 
cost function: 

e (x) = ^(A^(t)-^(x,i)) 2 / e 2 , (25) 

i,3 



where K^{t) denotes the experimental measurements of the 
decorrelation and the coefficients K, b j (x, t) are the theoretical 
decorrelations assuming that the defect is located at x. The 
typical fluctuations on the measured decorrelations are encap- 
sulated in the parameter e. 

To find the value of the scattering cross-section a, also un- 
known, we remark that e(x) is, as a function of a, a polyno- 
mial of degree two. There is therefore a minimum depending 
on x at 



CT ° P,(XJ ~c E,.vO«(x,*) a ■ 



(26) 



We reintroduce the value of <7 opt into the expression d25l > and 
get the optimized error function 



opt 



(*) = £ 






(27) 



which does not depend on a anymore. The most likely po- 
sition of the defect is the position xo of the minimum of e opt . 
The value of the cross-section is cr opt (xo) obtained from Equa- 
tion d26| i. 

To give an interpretation to the values of e(x), it is custom- 
ary to normalize it in the following way 



X™( x ) 



e(x) 
/ 



(28) 



where / = N — 4 is the number of degrees of freedom, 
since four model parameters -the cross-section and the carte- 
sian coordinates of the defect- are to be estimated. The quan- 
tity x^(x) has the following interpretations. If Xn( x ) ^" 1> 
it is very unlikely that the point xo is actually located at x. 
If Xn( x ) — 1 tne P om t x is a good candidate for x . If 
X^(xo) <C 1, there is a large area where x«( x ) < 1 which 
means that the inversion could not locate precisely the change 
because the value of e is too large. In other words the quality 
of measurements is too poor to give any satisfactory result. It 
is possible to use \\ ( x ) to obtain the probability density that 
the defect has appeared at the point x, which we define as: 



p(x) 



1 
C 



exp 



1 



;^x'( x ) 
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C 



exp 



e(x) 
2e 2 



(29) 



where C is a normalization constant such that J"p(x)dx = 1 
(see the appendix for a derivation of this formula). 




FIG. 7. Density of probability for the position of the moving scat- 
terer. 



at the center of the circle by adding a single scatterer with 
cross-section a. For each pair of receiver, we compute syn- 
thetic data through application of the formula (|2Qb . The Thou- 
less time td is defined as L 2 /D. As a measure of the pre- 
cision, a resolution length 6 is introduced, which we com- 
pute using the probability density function d29l ) as follows: 
8 2 = f(x — x ) 2 p(x)dx. 

In the vicinity of the change, we infer that the contribu- 
tions of the terms in e(x) are comparable and we deduce that 
S oc e. Thus, the precision with which the measurements are 
made directly influences the precision with which the change 
is located. We will not study the dependence of 8 with respect 
to e and we chose a value e = 0.01 throughout the numeri- 
cal study. Note that a uniform probability distribution corre- 
sponds to a complete absence of information concerning the 
location of the change, and gives the value 8 ~ L. The typical 
behaviour of the resolution 8 as a function of the number of 
source-receiver pairs is depicted in Figure|9] In the configura- 
tion described above, each pair gives a comparable contribu- 
tion to e(x) so that e(x) is approximately proportional to N. 
Therefore in the ideal case described in our example, we find 
that 8 oc N- 1 ' 2 . 

Note that the resolution cannot be made arbitrarily small by 
increasing N at will, because it is not possible to find an arbi- 
trary number of source-receiver pairs providing independent 
data. The value N entering into the scaling law 8 oc N~ 1 / 2 is 
the number of independent decorrelation measurements. 



B. Resolution versus number of source-receiver pairs 



C. Resolution versus coda time 



To investigate the resolution of the inversion technique de- 
pending on the parameters of the likelihood maximization, 
we use a numerical approach. We compute the best achiev- 
able resolution regardless of all experimental difficulties that 
potentially degrade the accuracy of the location. We use an 
ideal set-up made of one source and N receivers regularly 
distributed on a circle (see figure [HJ. We introduce a change 



The dependence of the resolution 8 with respect to coda 
time, is shown on Figure[lO]for a = 1, D = 1, c = 1, L = 10 
and e — 0.01. The resolution exhibits a minimum at a time of 
order i m ; n = tq ■ For a given source-receiver pair, the coda 
time is the time that has elapsed after the arrival of the bal- 
listic wave. Shortly after the ballistic arrival, the waves that 
reach the receiver have followed "snake-like" paths around 
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FIG. 8. Description of the numerical setup used for investigating 
the accuracy of the inversion technique. The example is shown with 
N = 5. The other parameters of the numerical simulations are : 
L = 10, D = 1, c = 1. The change xq is located at the center 
of the circle and is used in Section |TVB| to study the optimal spatial 
resolution of the inversion. The change located at point xi is used to 
study the robustness of the inversion technique against measurement 
errors on the determination of D in Section lTVEl 
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FIG. 9. Spatial resolution at the coda time t — td, as a function of 
the number of receivers, where td is the Thouless time. The typical 
setup for the numerical experiment is depicted in Figure[8] The dots 
correspond to the values of the resolution for N — 5, 10, 20, 40 
and 80. The double logarithmic scale provides clear evidence of the 
relation 8 ~ N~ 1/2 . 



the direct ray. In the early coda, the only signals sensitive to 
the change are those for which the change is located along 
the segment joining the source and the receiver. Later in the 
coda, the diffuse waves arriving at the receiver have explored a 
larger volume of the system. This qualitatively explains why 8 
decreases with the coda time t. At very late times, the formula 
(|20T > reveals that the decorrelation for each source receiver pair 
saturates, as the exponential factor tends to 1, The asymptotic 
spatial sensitivity to the change is algebraic only. After reach- 
ing a minimum, 8 increases because the variations of \n with 
respect to x decreases. The minimum for S is found approx- 
imately at time td, the Thouless time, after which the whole 
system has been explored by the diffuse waves and yet Q still 
exhibits large spatial variations. 
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FIG. 10. Spatial resolution obtained for the setup of Figure[8]for the 
values TV = 5, 10, 20, 40 and 80 and for coda times varying from 
2.10 -2 T£> to 10 3 td. The time scale is logarithmic. A minimum of 
the resolution is found at t ~ td- 
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FIG. 11. Optimal resolution 8 as a function of the cross-section of 
the change a obtained for the setup of Figure [8] Other parameters of 
the simulations are D — 1, c — 1, N — 10, e — 0.01. 



tion 8 decreases as a increases. Note that when a is very 
small, 8 goes to a value ~ L, meaning that it is not possi- 
ble to detect the change. When a ~ L 2 , the cross-section is 
equivalent to the area of the system, and locating a change has 
no physical significance in this limit. In Figure QT| we plot 
the variations of 8 at the optimal time t — td as a varies 
from 10 -4 I/ 2 to L 2 . The other parameters of the calculations 
are D = 1, c = 1, e = 0.01, N = 10. We observe that the 
spatial resolution 8 decreases by a factor 2 as the cross-section 
increases from 10~ 2 to 1. 



D. Resolution versus cross-section 

The scattering cross-section a of the change also influences 
the precision of the technique. We observe that the resolu- 



E. Sensitivity to the value of the diffusivity D 

Our inversion technique depends heavily upon our ability 
to estimate the diffusivity of the waves in the heterogeneous 
medium. Although the absorption time r does not enter into 




V. BOUNDARY CONDITIONS 
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FIG. 12. Effect of an estimation error on the value of the diffusivity 
D induced on the relocation of a target. Synthetic data were calcu- 
lated with D — 1 and inverted with modified values of the diffusivity 
D' ranging from 0.2 to 8. The other parameters of the simulation are 
a — 1, L — 10, N — 10, e = 0.01 and the change is located at 
the position xi (see Figure[8}. £ is the distance between xi and the 
point where \ 2 ( x ) is minimum. In the simulated configuration, the 
inversion technique remains accurate even if D' differs from D by a 
factor 2. 



The inversion technique presented in section HV Al relies on 
the knowledge of the function Pd, the diffusion kernel, which 
depends on the boundary conditions of the system. For sim- 
plicity, we studied the LOCADIFF technique in an infinite 
medium without taking into account the effect of boundaries, 
which may not be realistic in applications. An abundant liter- 
ature is dedicated to solving the diffusion equation in a wide 
range of situations [33]. In many cases of practical interest, 
sophisticated techniques are required to provide an exact so- 
lution or a numerical approximation up to a required accuracy. 
In the infinite medium, the decorrelation ( fTSI ) can be computed 
numerically. In the presence of boundaries, it is more difficult 
to compute the Green's function because translational invari- 
ance is lost. However, if the boundaries are fiat, it possible 
to construct the Green's function from the solution without 
boundaries using symmetry arguments. In the general case, 
one has to solve the diffusion equation for the geometry of the 
system, which is a problem for applied mathematics in itself. 

In the simple case of a single planar boundary, the solu- 
tion P ( f of the diffusion equation of the semi-infinite medium, 
can be deduced from PJ* 3 , using the technique of images: 



Pi(s, r, t) 



a(P d (s,r,t)+0P d {s',r,t)) 



(30) 



the final formula j20l ), let us remark that in practice D and 
t cannot be measured independently. The diffusivity D is 
the crucial physical parameter which enters into the formula 
for the intensity propagator Pd and controls the accuracy of 
the energy propagation model of the medium. It is therefore 
important to quantify the impact of errors in the diffusivity D 
on the accuracy of our method. Even if we use an incorrect 
value for the diffusivity, our inversion procedure still provides 
an answer for the position of the defect. The main issue is to 
quantify to what extent the inferred position differs from the 
exact location of the target. To address this point, we plot the 
spatial resolution and the absolute error of the inversion for a 
wide range of values of D on a specific example. 



We use the approach described in Section lTV Bl First, a syn- 
thetic data set is computed with a value D for the diffusivity. 
This synthetic data set is then inverted for the location of the 
target using a different diffusivity D'. The change is located 
at the position xi indicated on Figure [8] The other physical 
parameters L = 10, a = 1, N = 10, e = 0.01 and t have been 
adjusted to provide the smallest spatial resolution 6. We call £ 
the distance between the change located by the inversion and 
S is the resolution length. The results of the simulation are 
displayed in Figure Q~2] It is rather remarkable that an error 
on D as large as a factor of 2 yields a location of the change 
within one half of the resolution length. In this specific but 
realistic example, the inversion technique is therefore very ro- 
bust against errors on the determination of D. This constitutes 
a major advantage of our method. Based on these results, we 
infer that spatial variations of D within a factor of 2 will not 
affect the results dramatically. 



where s' is the image of s with respect to the boundary (see 
Figure [T3l > and j3 is a characteristic coefficient depending on 
the nature of the boundary condition. For instance if the 
boundary is absorbing, j3 = — 1 and if it is fully reflecting, 
we have (3 = 1. The normalization coefficient a is, in the case 
of constant diffusivity 



(31) 



where d,B, s is the distance from the source to the boundary. 
Note that a is undetermined in the case where the conditions 
ft = — 1 and d,B, s — are met simultaneously. The solution 





FIG. 13. Schematic representation of a boundary condition for Pd- 
The image of the source s is noted s' and the arbitrary point is r. The 
solid line is the solution in the infinite medium. 

to the diffusion equation in presence of the boundary can be 
plugged into the decorrelation expression ( TT8l > leading to four 
terms (figurePBli. Note that in case there are more boundaries, 
the image technique requires to take into account infinitely 
many images. Other techniques also lead to infinite series. 
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FIG. 14. In presence of a single straight boundary, the decorrelation 
function (TB) involves four terms, coming from the product of two 
formula of the form J30t . We use the images of the source s and the 
receiver r: Thanks to the symmetry of the diffusion equation, there 
is no need to introduce the image of the point x. 



VI. DISCUSSION 



In this section, we discuss issues related to the practical 
use of the LOCADIFF technique as well as possible improve- 
ments. We first note that if the interval between the records of 
hij and hL is large, the medium may also have experienced 
a global change, for instance a dilation due to a temperature 
change. In this case, the computation of the decorrelation may 
be refined by taking into account a global relative velocity 
change e, where e yields the maximum value of the correla- 
tion 



J hij((l + e)u) h'ij(u)du 

y/jhij^dujh^urdu'' 



(32) 



VII. CONCLUSION 

In this article, we have shown that it is possible to use the 
high sensitivity of diffuse waves to detect, characterize and 
locate a small change in a strongly scattering medium. Our 
technique uses the correlation of coda waveforms recorded 
before and after the change. Based on a maximum-likelihood 
approach, and a simple diffusion model, we demonstrate the 
possibility to retrieve the position of the change along with its 
scattering cross-section. We have also investigated the opti- 
mal values of the parameters that enter in the inversion proce- 
dure, based on a simple setup where sources and receivers are 
arranged on a circle surrounding the change. Three features 
have been identified: 1) We found that the resolution scales 
with the inverse square root of the number of sensors. 2) The 
technique provides the best results when the correlation win- 
dow is centered on the Thouless time of the system. 3) We 
demonstrated that the technique is not very sensitive to errors 
in the measurement of the diffusivity. 

Several aspects are still to be investigated. First, we have 
assumed that a single change occurs in the medium, an as- 
sumption which is probably too restrictive in some applica- 
tions. In a straightforward generalization of our technique to 
n changes, the dimension of the parameter space scales like 
An which in turn considerably increases the computation time. 
An alternative route for the inversion has to be found. Second, 
we have made the assumption of a point-like change. An ex- 
tended change may not necessarily be equivalent to a collec- 
tion of point-like changes. Again, an alternative approach to 
the inverse problem will be needed. We are currently investi- 
gating these two issues. 

Using 2D finite difference wave simulations, we have 
demonstrated that LOCADIFF efficiently locates a small 
change in a multiple scattering environment. In a seperate 
paper [34], experiments have also been conducted with ultra- 
sound in concrete. The change was a hole drilled in the sam- 
ple, and the LOCADIFF technique successfully retrieved its 
actual position. Other applications in geophysics and material 
sciences can be envisaged. 



where the integrals are performed along the whole record 11 

Another issue concerns the possible improvements on the 
inversion procedure. Under the form presented in this article, 
the LOCADIFF technique only uses a small time window in 
the signals. It would be of great interest to take into account 
several time windows in the coda. This would provide more 
independent data for the inversion procedure and may reduce 
the effect of noise. 

Finally, we point out that the kernel used in the inversion 
is computed from the solution to the diffusion equation. In 
some simple geometries, like the infinite medium, the solu- 
tion is analytic and simple to compute. If the shape of the 
medium is irregular, with possibly more complicated bound- 
ary conditions, the kernels can only be approximated numer- 
ically. Alternatively, our approach could benefit from recent 
developpements in the implementation of the radiative trans- 
fer equation. 
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Appendix A: Derivation using Bayesian inversion 

We shortly derive here the density of probability den- 
sity J29l using a Bayesian inference. In this calculation, we 
suppose that there is a change at an unknown position x. The 
values of the measurements i\~ 4 r ™ are accurate up to an error or- 
der e such that they are distributed around the numerical value 
Kij (x, i) according to a standard error function. 
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Each pair (i,j) provides an independent information. The 
Bayesian inversion consists in finding the probability den- 
sity of x knowing the values of Kfj, namely to compute 
p(x\{Kjj}). Let us call p n (x) the probability density for the 
position of the change when n source-receiver pairs have been 
taken into account. Before measurement, the probability of 
the location of the change is uniform in the whole medium, 
so we have po( x ) = h (V is the volume). Suppose we know 
p n - 1 (x) and let us compute the joint probability of x and K n 
using Bayes' formula. We use the two relations: 

p(x, K?\K?, . . . K™_ x ) = Pn (xMiC) (A2) 

p(x, K£\K?, . . . iC-i) = p(^|xK_i(x) (A3) 



Integrating ( IA21 > over x we can compute p(K™) as 

( j%„(x)dx) p(K?) = J p{x, K™\K™, . . . iC-i)dx 

(A4) 
The integral of p n (x) is equal to 1 so we conclude that, us- 
ing dA3l 



Pn(x) 



P (iC|xK-i(x) 
/ v P(^|xK-i(x)dx" 



(A5) 



Therefore we have a recurrence scheme yielding the distribu- 
tion of probability p N (x) : 



PnW 



;ntiP(^Tix)dx 



(A6) 



which gives Equation d29l i after replacing the probabilities 
with expression (IA1I) . 
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